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Modeling primary breakup: A three-dimensional 
Eulerian level set /vortex sheet method for 
two-phase interface dynamics 

By M. Herrmann 


1. Motivation and objectives 

Atomization processes play an important role in a wide variety of technical applications 
and natural phenomena, ranging from inkjet printers, gas turbines, direct injection IC- 
engines, and cryogenic rocket engines to ocean wave breaking and hydrothermal features. 
The atomization process of liquid jets and sheets is usually divided into two consecutive 
steps: the primary and the secondary breakup. During primary breakup, the liquid jet or 
sheet exhibits large scale coherent structures that interact with the gas-phase and break 
up into both large and small scale drops. During secondary breakup, these drops break 
up into ever smaller drops that finally may evaporate. 

Usually, the atomization process occurs in a turbulent environment, involving a wide 
range of time and length scales. Given today’s computational resources, the direct nu- 
merical simulation (DNS) of the turbulent breakup process as a whole, resolving all 
physical processes, is impossible, except for some very simple configurations. Instead, 
models describing the physics of the atomization process have to be employed. 

Various models have already been developed for the secondary breakup process. There, 
it can be assumed that the characteristic length scale £ of the drops is much smaller than 
the available grid resolution Aa: and that the liquid volume fraction in each grid cell 0; 
is small, see Fig. 1. Furthermore, assuming simple geometrical shapes of the individual 
drops, like spheres or ellipsoids, the interaction between these drops and the surrounding 
fluid can be taken into account. Statistical models describing the secondary breakup 
process in turbulent environments can thus be derived (O’Rourke 1981; O’Rourke & 
Amsden 1987; Reitz 1987; Reitz & Diwakar 1987; Tanner 1997). 

However, the above assumptions do not hold true for the primary breakup process. 
Here, the turbulent liquid fluid interacts with the surrounding turbulent gas-phase on 
scales larger than Ax, resulting in highly complex interface dynamics and individual grid 
cells that can be fully immersed in the liquid phase, compare Fig. 1. An explicit treatment 
of the phase interface and its dynamics is therefore required. To this end, we propose to 
follow in essence a Large Eddy Simulation (LES) type approach: all interface dynamics 
and physical processes occurring on scales larger than the available grid resolution Aa: 
shall be fully resolved and all dynamics and processes occurring on subgrid scales shall 
be modeled. The resulting approach is called Large Surface Structure (LSS) model. 

In order to develop such a LSS model for the turbulent primary breakup process, 
one potential approach is to start off from a fully resolved description of the interface 
dynamics using the Navier-Stokes equations and include an additional source term in the 
momentum equation due to surface tension forces (Brackbill et al. 1992). In order to track 
the location, motion, and topology of the phase interface, the Navier-Stokes equations are 
then coupled to one of various possible tracking methods, for example marker particles 
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Figure 1. Breakup of a liquid jet. 


(Brackbill et al. 1988; Rider & Kothe 1995; Unverdi & Tryggvason 1992), the Volume-of- 
Fluid method (Noh & Woodward 1976; Kothe & Rider 1994; Gueyffier et al. 1999), or the 
level set method (Osher & Sethian 1988; Sussman et al. 1994, 1998). Then, introducing 
ensemble averaging or spatial filtering results in unclosed terms that require modeling 
(Brocchini & Peregrine 2001a, 6). Unfortunately, the derivation of such closure models is 
not straightforward and, hence, has not been achieved yet. This is in part due to the fact 
that, with the exception of the surface tension term, all other physical processes occurring 
at the phase interface itself, like for example stretching, are not described by explicit 
source terms. Instead, they are hidden within the interdependence between the Navier- 
Stokes equations and the respective interface tracking equation. Thus, a formulation 
containing the source terms explicitly could greatly facilitate any attempt to derive the 
appropriate closure models. 

To this end, a novel three-dimensional Eulerian level set/vortex sheet method is pro- 
posed. Its advantage is the fact that it contains explicit source terms for each individual 
physical process that occurs at the phase interface. It thus constitutes a promising frame- 
work for the derivation of the LSS subgrid closure models. 

This paper is divided into four parts. First, the level set/vortex sheet method for three- 
dimensional two-phase interface dynamics is presented. Second, the LSS model for the 
primary breakup of turbulent liquid jets and sheets is outlined and all terms requiring 
subgrid modeling are identified. Then, preliminary three-dimensional results of the level 
set/vortex sheet method are presented and discussed. Finally, conclusions are drawn and 
an outlook to future work is given. 

2. The level set/vortex sheet method 

The aim of the level set/vortex sheet method is to describe the dynamics of the phase 
interface T between two inviscid, incompressible fluids 1 and 2, as shown in Fig. 2. In this 
case, the velocity on either side i of the interface T is determined by the incompressible 
Euler equations, given here in dimensionless form, 



V • Ui = 0 , 


1 


Vp , 


(2.1) 

( 2 . 2 ) 


Pi 
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subjected to the boundary conditions at the interface T, 


[(«i ^ “ 2 ) • n] 
[n x (u 2 - «i)] 

\P2 - Pi] 


= 0 
r 



r We* 


(2.3) 

(2.4) 

(2.5) 


and at infinity, 


lim Ui = ±«oo . (2-6) 

y — >-=hoc 

Here, n is the interface normal vector, r] is the vortex sheet strength, and k is the local 
curvature of T. The Weber number is defined as 


We = PrefMref/^ref , (2.7) 

where S is the surface tension coefficient and p re f, u re f, and L re f are the reference den- 
sity, velocity, and length, respectively. An interface subjected to the above boundary 
conditions is called a vortex sheet (Saffman & Baker 1979). 

The partial differential equation describing the evolution of the vortex sheet strength r) 
can be derived by combining the Euler equations, Eqs. (2.1) and (2.2), with the boundary 
conditions at the interface, Eqs. (2.3)-(2.5), resulting in (Pozrikidis 2000) 

dri 

+ u ■ X7rj = —n x [( r) x n) • Vit] + n [(Vm • n) ■ rj\ 

2(a + lj 

H — (n x Vk) + 2 An x a . (2.8) 

We 

Here, A = (pi — p 2 )/(pi + p 2 ) is the Atwood number and a is the average acceleration 
of fluid 1 and fluid 2 at the interface. The major advantage of Eq. (2.8), as compared to 
a formulation based on the Euler equations, is the fact that Eq. (2.8) contains explicit 
local individual source terms on the right-hand side describing the physical processes at 
the interface. These are, from left to right, two stretching terms, a surface tension term 
T a , and a density difference term. 

In addition to the evolution of the local vortex sheet strength, Eq. (2.8), the location 
and motion of the phase interface itself has to be known. To this end, vortex sheets are 
typically solved by a boundary integral method within a Lagrangian framework where 
the phase interface is tracked by marker particles (Baker et al. 1982; Pullin 1982; Hou 
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et al. 1997, 2001; Rangel & Sirignano 1988). Marker particles allow for highly accurate 
tracking of the phase interface motion in a DNS. However, the introduction of ensemble 
averaging and spatial filtering of the interface topology is not straightforward and hence 
a strategy for the derivation of appropriate LSS subgrid closure models is not directly 
apparent. 

Level sets, on the other hand, have been successfully applied to the derivation of closure 
models in the field of premixed turbulent combustion (Peters 1999, 2000). Thus, instead 
of using marker particles to describe the location and motion of the phase interface, here, 
the interface is represented by an iso-surface of the level set scalar field G(x, t), as shown 
in Fig. 2. Setting 

G(x, t) | r = Go = const , (2.9) 

G(x,t ) > Go in fluid 1, and G(x,t ) < Go in fluid 2, an evolution equation for the scalar 
G can be derived by simply differentiating Eq. (2.9) with respect to time, 

DG 

— + «'VG = 0. (2.10) 

This equation is called the level set equation (Osher & Sethian 1988). Using the level set 
scalar, geometrical properties of the interface, like its normal vector and curvature, can 
be easily expressed as 

VG ~ (2.11) 


n = 


|VG| ’ 


k = V • n . 


Strictly speaking, Eqs. (2.8) and (2.10) are valid only at the location of the inter- 
face itself. However, to facilitate the numerical solution of both equations in the whole 
computational domain, T) is set constant in the interface normal direction, 


V»7 • VG = 0 , 

and G is chosen to be a distance function away from the interface, 


|VG| 


GjtG o 


= 1 . 


(2.12) 


(2.13) 


Equations (2.8) and (2.10) are coupled by the self-induced velocity u of the vortex 
sheet. To calculate u, the vector potential ip is introduced, 


Axp = u) . 


(2.14) 


Here, the vorticity vector u; is calculated following a vortex-in-cell type approach (Chris- 
tiansen 1973; Cottet & Koumoutsakos 2000) 

oj{x) = I t){x')5{x — x')5 {G{x') — Go)\S7G{x')\dx' , (2.15) 

Jv 


where 8 is the delta-function. Then, u can be calculated from 


u(x) = [ 8(x — x') (V x xp) dx' . 

Jv 


(2.16) 


In summary, Eqs. (2.8), (2.10), and (2.14) - (2.16) describe the three-dimensional two- 
phase interface dynamics and constitute the level set/vortex sheet method. 

2.1. Numerical methods 

Numerically, Eqs. (2.8) and (2.10) are solved in a narrow band (Peng et al. 1999) by 
a third-order WENO scheme (Jiang & Peng 2000) using a third-order TVD Runge- 
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Figure 3. LSS interface definitions. 


Kutta time discretization (Shu & Osher 1989). The redistribution of r] (2.12) is solved 
by a Fast Marching Method (Sethian 1996; Adalsteinsson & Sethian 1999), whereas the 
reinitialization of G (2.13) is solved by an iterative procedure (Sussman et al. 1994; 
Peng et al. 1999). The interested reader is referred to Herrmann (2002) for a detailed 
description of the numerical methods employed in the level set/vortex sheet method. 


3. The LSS model for turbulent primary breakup 

The basic idea of the LSS model is to split the treatment of the primary breakup 
process into two parts. All phase interface dynamics occurring on scales larger than 
the local grid size are explicitly resolved and tracked by a level set approach, whereas 
interface dynamics occurring on subgrid scales are described by an appropriate subgrid 
model. Furthermore, the LSS subgrid model has to separate out all broken off subgrid 
scale liquid drops and transfer them to a secondary breakup model. 

The level set equation describing the interface location and motion on the resolved 
scales can be derived by first introducing appropriate interface based filters (Oberlack 
et al. 2001) into the level set equation (2.10), see Fig. 3, 

dG ^ 

— + mAG = 0. (3.1) 

Here, ~ denotes quantities on the resolved (filter) scale. Furthermore, the mass transfer 
rate rh p into the secondary breakup model has to be taken into account, 

m P = Pis p A Go = ^7rpiJ^ J P(D)D 3 dD , (3.2) 

where s p is the subgrid primary breakup velocity, Aq 0 is the local surface area of the 
resolved interface, and P(D) is the droplet diameter number distribution, 

P{D) = N , (3.3) 

where N is the total number of drops. Then, the resolved scale level set equation reads 



dG 

qj + (u + s p n) ■ VG = 0 , 


(3.4) 
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where n is the normal vector of the resolved scale interface 


VG 

n = — — . 
|VG| 


(3.5) 


To describe the phase interface dynamics on the resolved scale, their effect on the flow 
field has to be taken into account by an additional source term T in the momentum 
equation, 

T = T,k5(G — Go)n + T SGS . (3.6) 

Here, the first term on the right-hand side describes the effect of surface tension forces 
due to the local curvature k, of the resolved scale interface, whereas the second term, 
T sgs 5 accounts for the effect of the subgrid scale surface tension forces on the resolved 
scale flow field. 

Thus, the yet unclosed subgrid terms of the LSS model requiring modeling are the 
subgrid primary breakup velocity s p , the droplet diameter number distribution P(D), 
and the subgrid scale surface tension effect Tsgs- As previously indicated, these subgrid 
terms are to be derived from the level set/vortex sheet method. Performing DNS of 
the primary breakup of liquid surfaces and sheets in turbulent environments will help 
to identify characteristic regimes of the turbulent primary breakup and their dominant 
physical processes. These can then be quantified using the explicit source terms in the 
77 -equation (2.8), thus providing guidelines for the derivation of appropriate LSS subgrid 
models. 


4. Results 

In order to both validate the three-dimensional level set/ vortex sheet method and to 
demonstrate its ability to perform DNS of the primary breakup process, the results of two 
different cases are presented. First, the calculated oscillation periods of liquid columns 
and spheres are compared to theoretical results. Then, the breakups of a randomly per- 
turbed liquid surface and sheet are presented. 

4.1. Oscillating columns and spheres 

To validate the proposed level set/vortex sheet method, the calculated oscillation periods 
T of liquid columns and spheres of mean radius R = 0.25, center x c = (0.5, 0.5, 0.5), 
amplitude e = 0.051?, and Atwood number A — 0 are compared to theoretical results 
(Lamb 1945). The initial vortex sheet strength in both cases is set to 

r/(x, t = 0) = 0 . (4.1) 

All calculations are performed in a unit sized box resolved by an equidistant cartesian 
grid of 128 x 128 and 128 x 128 x 128 nodes, respectively. 

Figure 4 shows the distribution of the surface tension term T a of the 77 -equation (2.8) 
in the x-, y-, and 3-direction, 

Tct = (n x Vk) ’ (4 - 2) 

for the oscillating sphere of mode number n = 5 and Weber number We = 10 calculated 
at t = 0. As the shape of the sphere indicates, T a in the ^-direction is a factor of roughly 
four higher than T a in the other two directions, leading to the predominant oscillation 
in the y- 3 -plane. 
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Figure 4. Distribution of the surface tension term T a in the x-direction (top), y-direction 
(left), and ^-direction (right) for the mode n = 5 oscillating sphere at t = 0 and We = 10. 


Figure 5 depicts the comparison of the oscillation period for the oscillating columns 
on the left-hand side and the oscillating spheres on the right-hand side for two different 
Weber numbers. As can be clearly seen, agreement between simulation and theory is very 
good. 

4.2. Liquid surface and sheet breakup 

To demonstrate the capability of the proposed level set / vortex sheet method to simulate 
the primary breakup process, the temporal evolution of both a randomly perturbed 
liquid surface and sheet are simulated. In the case of the liquid surface, the on average 
flat interface located at z = 0 is perturbed in the ^-direction by a Fourier series of 64 
sinusoidal waves in both the x- and y-direction with random amplitude 0 < e < 0.01 and 
random phase shift. In the case of the liquid sheet, the two on average flat interfaces are 
located at z = —B / 2 and z = +B/2 and are again perturbed by two Fourier series of 64 
sinusoidal waves. The thickness of the liquid sheet is set to B = 0.1. 

The initial vortex sheet strength for the liquid surface is set to 

r/(x,t = 0) = (—1,0,0) (4.3) 


and to 


T]{x,t = 0) 


(-1,0,0) : z > 0 

( 1,0,0) : z<0, 


(4.4) 
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Figure 5. Oscillation period T of liquid columns (left) and spheres (right) as a function of mode 
number n for varying Weber numbers We. Lines denote theoretical and symbols computational 
results. 


in the liquid sheet case. Both, the surface and the sheet simulation were performed in a 
x- and y-direction periodic box of size [0, 1] x [0, 1] x [—1,1] resolved by a cartesian grid 
of 64 x 64 x 128 equidistant nodes. In both simulations, the Atwood number is A = 0. 
The Weber number in the surface simulation is We = 500 and the Weber number in the 
sheet simulation based on the sheet thickness is Wes = 100. 

As depicted in Fig. 6, the surface shows an initial growth of two-dimensional Kelvin- 
Helmholtz instabilities (t = 1). These continue to grow ( t = 3) and form three-dimensional 
structures ( t = 5) resulting in elongated fingers ( t = 6.5) that finally initiate breakup 
( t = 8 . 0 ). 

The liquid sheet, depicted in Fig. 7, also exhibits the initial formation of two-dimensional 
Kelvin-Helmholtz instabilities (t = 1) that continue to grow (t = 3) until the liquid film 
gets too thin and ruptures (t = 5). Individual fingers are formed that extend mostly in 
the transverse direction (t = 8) and continue to break up into individual drops of varying 
sizes (t = 12). 


5. Conclusions and future work 

A Eulerian level set/ vortex sheet method has been presented that allows for the three- 
dimensional calculation of the phase interface dynamics between two inviscid and in- 
compressible fluids. Results obtained with the proposed method for oscillating columns 
and spheres show very good agreement with theoretical predictions. Furthermore, the 
applicability of the method to the primary breakup process has been demonstrated by 
simulations of the breakup of both a liquid surface and a liquid sheet. 

In addition, the LSS model for turbulent primary breakup has been outlined, and all 
terms requiring subgrid modeling have been identified. The proposed level set/vortex 
sheet method has the advantage that it allows for the detailed study of each individual 
physical process occurring at the phase interface. It thus provides a promising framework 
for the derivation of the LSS subgrid models. 

Future work will focus on including the effect of non-zero Atwood numbers and on 
coupling of the level set/vortex sheet method to an outside turbulent flow field. Also, the 
level set/vortex sheet method will be parallelized making use of the new domain decom- 
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Figure 6. Temporal evolution of the three-dimensional liquid surface breakup, A = 0, 

We = 500. 


position parallelization of the Fast Marching Method presented in Herrmann (2003). This 
will allow for efficient DNS of the primary breakup process to help identify the different 
regimes of turbulent primary breakup and their dominant physical processes, facilitating 
the derivation of the LSS subgrid models. Finally, combining the LSS model to spray 
models describing the secondary breakup will allow for the first LES of the turbulent 
atomization process as a whole. 
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11 11 

Figure 7. Temporal evolution of the three-dimensional liquid sheet breakup, A = 0, 

Wes = 100. 
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